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20.  ABSTRACT  (Continued) 

capabilities  of  the  NRL  reverberation  model  are:  the  treatment  of  both  monostatic  and  bistatic 
geometries;  the  inclusion  of  all  multipath  contributions  in  the  reverberation  calculation;  range- 
dependent  propagation  effects  (presently  for  monostatic  geometries);  angle-dependent  surface 
and  bottom  scattering  strength  functions;  frequency  spreading  of  reverberation  returns  induced 
by  forward  scattering  at  the  sea  surface;  and  the  doppler-shifting  of  surface  backscattered  energy. 
The  frequency  spreading  and  doppler-shifting  estimates  are  calculated  by  supplemental  programs 
which  are  not  presently  included  in  the  core  reverberation  program  to  be  described  in  this  report. 

The  reverberation  model  comprises  a  sequence  of  computer  programs.  This  allows  the 
evaluation  of  intermediate  outputs,  and  facilitates  the  model’s  application  to  parameter  studies 
by  delaying  parameter  inputs  until  the  calculation  stage  in  which  they  are  needed  is  reached.  In 
the  first  stage,  sound  speed  profiles  and  bottom  elevations  are  input  as  a  function  of  range.  Nor¬ 
mal  dimensioning  permits  for  up  to  SO  sound  speed  profiles  and  SO  bottom  elevations  but  is 
easily  expanded  to  cover  more  refined  environmental  specification.  Next,  a  ray-acoustic  propaga¬ 
tion  model  is  used  to  trace  out  rays  from  both  source  and  receiver,  keeping  track  of  various 
parameters  at  each  boundary  encounter  for  each  ray  path.  Signal  intensity  is  modified  by  bottom 
losses,  geometrical  spreading,  and  volume  absorption  which  varies  as  a  function  of  frequency. 
The  number  of  acoustic  rays  traced  is  a  variable  and  is  selected  to  ensure  adequate  sampling  of 
the  propagation  field.  Then,  the  reverberation  contribution  from  each  elemental  scattering  area 
on  the  boundary  is  obtained  by  calculating  the  time-dependent  contributions  due  to  all  possible 
reverberant  paths  from  source  to  scattering  area  to  receiver.  That  is,  the  calculation  utilizes  all 
m  outgoing  paths  combined  with  n  backscattered  return  paths.  Each  contribution  is  weighted  by 
source  and  receiver  beam  patterns,  the  angle-dependent  surface  or  bottom  backscattering 
strength  values,  and  the  duration  and  shape  of  the  source  pulse.  Finally,  by  integrating  over  the 
entire  scattering  surface,  the  total  boundary  reverberation  is  determined.  The  main  outputs  of 
the  reverberation  model  are  average  reverberation  vs  time  plots  and  average  reverberation  den¬ 
sity  as  a  function  of  vertical  angle  at  selected  times.  When  frequency  spreading  and  doppler- 
shifting  calculations  are  included,  it  is  also  possible  to  plot  reverberation  spectral-time  histories. 

2  py ? o tp d  «i4 

"Tills  leport  pieseuts  the  theoretical  foundations  of  the  reverberation  model,  the  numerical 
implementation  of  this  theory,  and  a  detailed  description  of  actual  program  execution.  A  com¬ 
plete  description  of  all  inputs  and  outputs  is  included.  Finally,  examples  illustrating  the  use  of 
the  sequence  of  programs  are  provided. 
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NRL  REVERBERATION  MODEL:  A  COMPUTER  PROGRAM 
FOR  THE  PREDICTION  AND  ANALYSIS  OF  MEDIUM-  TO 
LONG-RANGE  BOUNDARY  REVERBERATION 


INTRODUCTION  ^ '  A  « -  -  ~ 

A  sequence  of  computer  programs^  to  predict  long-range,  Lotf"  frequency  monostatic  or  bistatic 
boundary  reverberation  from  either  the  ocean  surface  or  bottom Ihas  been  developed.  The  formulation 
requires  sets  of  acoustic  rays  to  be  traced  out  from  both  the  source  and  receiver,  keeping  track  of  vari¬ 
ous  parameters  at  erch  boundary  encounter.  The  reverberation  contribution  from  each  elemental 
scattering  area  on  the  boundary  is  obtained  by  calculating  the  time-dependent  contributions  due  to  all 
possible  reverberant  paths  from  source  to  scattering  area  to  receiver.  That  is,  the  calculation  utilizes  all 
outgoing  paths  combined  with  all  backscattered  incoming  paths.  Each  contribution  is  weighted  by  the 
source  and  receiver  beampatterns,  the  ang’e-dependent  surface  or  bottom  backscattering  strength 
values,  and  the  duration  and  shape  of  the  source  pulse.  Finally,  by  integrating  over  the  entire  scatter¬ 
ing  surface,  the  total  boundary  reverberation  is  determined.  Some  of  the  more  significant  features  of 
this  computer  model  are:  its  ability  to  treat  both  monostatic  and  bistatic  geometries;  the  capability  of 
modeling  the  ocean  environment  as  range  dependent  when  the  source  and  receiver  are  not  separated 
horizontally;  the  ability  to  input  source  and  receiver  positions  and  beam  patterns;  the  use  of  angle- 
dependent  scattering  coefficients  at  the  surface  and  bottom;  and  the  ability  to  select  the  duration  and 
shape  of  the  source  pulse.  It  is  also  possible  to  include  frequency  spreading  from  forward  surface 
scattering  and  doppler-shifted  surface  backscattering  calculations.  At  present,  the  model  is  capable  of 
handling  a  true  bistatic  geometry  for  a  range-independent  environment  only.  The  main  outputs  of  the 
model  are  plots  of  reverberation  vs  time,  and  reverberation  as  a  function  of  vertical  angle  at  selected 
times.  Table  1  summarizes  the  input  parameters  for  both  monostotic  and  bistatic  reverberation  calcula¬ 
tions. 


The  theoretical  model  is  discussed  in  the  first  section.  Computer  implementation  is  then 
described.  The  development  is  similar  to  that  given  in  a  previous  report  [11.  The  remainder  of  the 
report  describes  in  detail  the  means  of  actually  executing  the  computer  program.  A  complete  descrip¬ 
tion  of  all  inputs  and  outputs  is  included.  Finally,  examples  illustrating  the  use  of  the  sequence  of  pro¬ 
grams  are  provided. 


Table  1  —  Reverberation  Model  Parameters 


MONOSTATIC 

BISTATIC 

Source/Receiver  Positions 

Source/Receiver  Patterns 

Range  Dependent  Sound  Speed 

Range  Dependent  Bottom  Profile 

Angle  Dependent  Bottom  Loss  Tables 

Angle  Dependent  Surface  Scattering  Strength 
Angle  Dependent  Bottom  Scattering  Strength 
Pulse  Length,  Selectable 

Same 

Same 

Constant  Profile 
Flat  Bottom 
Same 

Same 

Same 

Same 
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THEORETICAL  MODEL 

To  give  a  qualitative  description  of  the  way  in  which  boundary  reverberation  arises,  suppose  that  a 
point  source  S  radiates  a  time-dependent  signal  into  the  three-dimensional  ocean.  The  acoustic  energy 
is  propagated  away  from  S  along  the  ray  paths.  When  a  ray  encounters  an  ocean  boundary,  it  is 
assumed  to  continue  to  propagate  in  the  direction  of  specular  reflection,  with  a  small  amount  of  the 
incident  radiation  being  reradiated  or  scattered  according  to  some  scattering  law.  The  scattering  is 
assumed  to  arise  from  the  excitation  of  small  scattering  elements  in  the  boundary,  each  of  which  acts  as 
a  weak  secondary  point  source  of  radiation.  Some  of  the  scattered  rays  will  arrive  at  the  receiving  point 
R,  and  the  sum  of  these  returns,  as  a  function  of  time,  is  the  boundary  reverberation. 

Because  the  frequency  of  the  reverberation  arising  from  surface  scattering  will  usually  be  doppler 
shifted,  and  that  from  bottom  scattering  will  not,  we  can  calculate  these  two  types  of  reverberation 
separately.  Thus  for  either  type  of  reverberation  it  will  be  necessary  to  trace  a  set  of  acoustic  rays  out 
from  both  5  and  R,  keeping  account  of  travel  times,  ray  intensities,  and  other  parameters  at  each 
boundary  encounter.  (If  the  receiving  point  R  is  coincident  with  the  source  point  S,  only  one  ray  trace 
is  necessary.)  To  develop  a  model  of  the  reverberation  process  from  which  an  expression  for  the  mean 
reverberation  decay  can  be  deduced,  three  specific  assumptions  concerning  the  scattering  process  are 
made: 

1.  A  small  element  of  area  dA  on  the  boundary,  when  excited  by  a  wave  that  is  nearly  plane, 
acts  as  a  secondary  point  source  for  exactly  the  duration  it  is  excited  by  the  incident  wave. 
Furthermore,  the  intensity  of  this  secondary  point  source  is  proportional  to  the  incident 
intensity  and  to  dA,  the  constant  of  proportionality  being  the  scattering  strength  <r. 

2.  For  purposes  of  calculating  a  mean  reverberation  envelope,  interference  effects  associated 
with  differences  in  acoustic  phase  can  be  ignored.  Thus  each  ray  arrival  at  a  scattering  ele¬ 
ment  produces  a  set  of  scattered  rays,  each  having,  in  effect,  a  random  phase  shift  relative  to 
the  incident  ray.  This  leads  to  the  determination  of  a  mean  reverberation  envelope  which  is 
representative  of  ensemble-averaged  reverberation  returns.  An  envelope  based  on  coherent 
summation  of  scattered  paths  would  be  more  indicative  of  a  single  sample  within  an  ensem¬ 
ble. 

3.  The  scattering  layer  at  the  ocean  surface  can  be  approximated  by  a  horizontal  plane  and 
appropriate  scattering  coefficient,  and  similarly  at  the  bottom  the  scattering  surface  has 
approximately  the  same  topography  as  the  bottom,  with  its  appropriate  scattering  coefficient. 

The  basic  procedure  for  calculating  reverberation  when  employing  the  concept  of  scattering 
strength  is  straightforward.  For  example,  at  an  element  of  the  ocean  surface  from  which  surface  rever¬ 
beration  is  to  be  calculated,  there  is  associated  a  set  of  reverberant  paths,  each  making  an  elemental 
contribution  to  the  reverberation.  Each  elemental  contribution  is  composed  of  the  transmission  losses 
associated  with  its  reverberant  path  from  the  source  to  the  scattering  point  and  from  the  scattering 
point  to  the  receiver,  appropriately  weighted  by  the  source  and  receiver  beam  patterns,  the  scattering 
strength,  and  the  area  of  the  element.  The  expected  or  average  value  of  the  instantaneous  reverbera¬ 
tion  is  then  the  sum  of  all  those  elemental  contributions  active  at  that  instant. 

A  general  expression  for  calculating  boundary  reverberation  can  be  formulated  based  on  the 
assumptions  outlined  above.  We  will  first  derive  the  expression  for  monostatic  boundary  reverberation. 
The  extension  to  bistatic  reverberation  will  then  be  made. 

Monostatic  Boundary  Reverberation 

Figure  1  shows  a  source-receiver  S-R  that  ensonifies  a  scattering  element  dA  via  a  propagation 
path  i  and  receives  reverberation  from  dA  via,  in  this  case,  a  different  return  path  j.  (The  figure  de¬ 
picts  surface  reverberation.  For  bottom  reverberation  the  coordinates  of  dA  would  be  (X,  Y,  Zh)  where 
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Fig.  1  —  Source-receiver  ensonifying  a  scattering  element  dA  via  propagation 
path  /  and  receiving  reverberation  via  return  path  j 

Zb  is  the  ocean  depth.)  It  is  supposed  that  S  emits  a  time-varying  signal  of  intensity  /i(f),  weighted  by 
a  beam  pattern  B.  This  intensity  is  diminished  by  the  transmission  loss  L,  along  path  /  from  S  to  dA. 
The  incident  intensity  at  dA  is  then  weighted  by  the  scattering  strength  a.  The  reverberant  return  is 
further  diminished  by  the  transmission  loss  Lj  along  path  j  from  dA  to  R  At  R  it  is  weighted  by  the 
receiver  beam  pattern  B.  Further,  the  reverberation  return  is  delayed  in  time  by  the  two-way  travel 
time  Tj  +  7}.  That  is,  the  unweighted  initial  source  intensity  which  gives  rise  to  the  reverberation  con¬ 
tribution  at  time  t  due  only  to  outward  path  i  and  return  path  j  is  I\(t  —  T,—  Tj).  By  summing  over  all 
the  multipath  combinations  between  S—R  and  dA,  the  total  contribution  from  dA  is  obtained.  These 
intensities  will  be  summed  with  regard  to  time  but  without  regard  to  phase,  owing  to  the  second 
assumption.  Finally  by  integrating  over  the  area  A  of  scattering  points,  the  total  reverberation  R  (r )  at 
time  i  is  given  by 

n-r/x,  »>  turnin'* (1) 

where  m(X,Y)  indexes  the  set  of  rays  between  S—R  and  dA.  Note  that  the  scattering  strength 
cr (0/,  4>it  9'j,  <f>j)  in  Eq.  (1)  depends  on  the  inclination  and  azimuth  angles  of  incidence  and  scatter; 
that  is,  at  this  level  of  generality,  a  full  three-dimensional  scattering  pattern  would  be  employed. 

It  is  not  possible  to  numerically  evaluate  Eq.  (1)  directly  due  to  the  lack  of  a  three-dimensional 
propagation  model.  Therefore  the  monostatic  reverberation  model  was  reduced  to  two  dimensions  by 
assuming 

1.  cylindrical  symmetry  of  the  ocean  environment  about  the  Z-axis  (with  the  cylindrical 
coordinates  being  Gc ,z,<f>)) 

2.  directions  of  reflected  and  scattered  energy  are  confined  to  the  vertical  x  -  z  plane  con¬ 
taining  the  incident  ray  path. 

It  is  also  assumed  that  the  beam  patterns  of  source  and  receiver  are  both  separable,  and  we  write 

5(0,  <f>)  -  a(d>)b(0),  (2) 


B(9,  0)  -  a(<f>)  b(9). 


~ajr 
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Let 


<*>=“(*  a(4>)a(4>)  dj> 

*'0 

be  the  equivalent  azimuthal  beamwidth  of  S-R.  Then  Eq.  (1)  becomes 

*(,)  -  *  £  ,J„  ,i('  -  -  rA>)  & 


(4) 


(5) 


Next,  we  take  the  signal  emitted  by  S  to  be  a  wave  of  constant  intensity  I\  and  duration  D  that  is 
turned  on  at  5  at  t  —  0.  Then  only  those  reverberating  elements  at  ranges  x  satisfying 


0  <  t  -  T,(x)  -  Tjix)  <  D 


(6) 


are  active  in  producing  reverberation  at  R  at  time  t.  The  double  inequality  of  Eq.  (6)  is  equivalent  to 

t  -  D  <  7;(x)  +  Tj(x)  <  t.  (7) 


Therefore  the  set  XyU)  of  reverberators  active  via  the  path  pair  (i,  j)  at  time  /  is  defined  by 

X,j(t)  =  xt  -  D  <  Ti(x)  +  Tj{x)  <  t 

The  characteristic  function  C\  of  a  set  \  may  be  defined  by 


C,(x) 


1  x 
o  X9x 


The  reverberation  R  (t)  at  R  can  then  be  written  as 


*<,)  -  '■*  A,  Ju,  cv *  “*■ 


(8) 


(9) 


(10) 


It  has  been  shown  [2]  that  this  formulation  reduces  to  the  standard  result  [3]  when  the  sound 
speed  c  and  scattering  coefficient  a-  are  constant  and  the  source-receiver  is  placed  at  a  depth  z0  below 
the  surface  in  a  bottomless  ocean.  That  result,  for  times  t  >>  D  +  2 (z0/c),  is 

R(t)  ~  M  .  (11) 


cd| 

ct 

2 

2 

In  terms  of  the  path  length  /  -  ct/ 2, 


/?(/)  ~  «7- 


cD , 

—  I  I* 


(12) 


the  quantity  in  the  brackets  being  the  reverberating  area.  It  is  straightforward  to  extend  this  result  to 
the  case  of  multiple  reflections  in  a  channel  with  a  bottom  at  a  finite  depth. 

Equation  (10)  is  a  generalization  of  the  isospeed  result  within  the  limitations  of  the  assumptions 
made  in  deriving  it.  The  generalization  is  twofold.  First,  Eq.  (10)  can  accommodate  a  wide  class  of 
sound-speed  fields,  including  horizontal  as  well  as  vertical  gradients,  provided  the  assumption  of 
cylindrical  symmetry  is  reasonably  well  met  over  the  significant  azimuthal  directions  defined  by  the 
source-receiver  beam  patterns.  Thus  refractive  effects  are  taken  into  account.  Second,  reverberation 
arising  from  different  outgoing  and  incoming  paths  is  included  in  addition  to  reverberation  arising  only 
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from  coincident  outgoing  and  incoming  paths.  At  shorter  ranges  the  "noncoincident"  contributions  are 
probably  not  appreciable,  but  at  longer  ranges  they  may  be  significant  even  with  vertically  directive 
beam  patterns. 

Bistatic  Boundary  Reverberation 

Equation  (1)  may  be  generalized  to  the  case  where  the  source  and  receiver  are  separated  in  range 
(the  bistatic  case).  The  times  at  which  an  elemental  contribution  of  the  boundary  to  bistatic  reverbera¬ 
tion  is  active  are  determined  by  the  travel  times  along  the  propagation  and  return  paths  and  the  dura¬ 
tion  of  the  transmitted  pulse.  Again  we  take  the  emitted  signal  to  be  a  wave  of  constant  intensity  /| 
and  duration  D  that  is  turned  on  at  S  at  t  -  0.  Then,  a  scattering  point  is  actively  contributing  at  time  t 
via  a  propagation  path  i  and  a  return  path  j  if 

t  -  D  <  T,(X,Y)  +  Tj(X.Y)  <  i,  (13) 

where  Tt  and  7}  are  the  one  way  travel  times  along  the  paths  /'  and  j.  The  set  of  points  active  at  time  / 
via  the  path  pair  (i.j)  is  then: 

A0U)  =  (( X,Y)\t-D  <  Ti(X.Y)  +  Tj(X,Y)  <  f).  (14) 


The  reverberation  R(t)  from  the  boundary  area  A  can  then  be  written: 


*  w  -  h  SSA  £  I  CA  (. X,  Y)  -{°{xY)llxY)~  a(9!>  +1’  V’  ^  (15) 

*  Itm(X.Y)  j€n(X,Y)  V  LjKA,  I  )  Lj\A,  I ) 


where  the  characteristic  function  CA  (X,Y)  is  1  if  (X,  K)  €  A0(t)  and  0  otherwise,  B  is  the  source 
beam  pattern  applied  to  path  /,  L,  is  the  accumulated  loss  along  path  /'  to  the  scattering  point,  Lj  is  the 
accumulated  loss  along  path  j  from  the  scattering  point  to  the  receiver,  B  is  the  receiver  beam  pattern 
applied  to  path  j,  and  <r  is  a  fully  three-dimensional  scattering  strength.  Note  that,  in  Eq.  (15), 
m(X,Y)  indexes  the  set  of  rays  between  the  source  5  and  dA,  while  n(X,Y)  indexes  the  set  of  rays 
between  the  receiver  R  and  dA.  For  a  bistatic  reverberation  calculation,  these  two  index  sets  are  not 
necessarily  the  same.  This  reflects  the  fact  that  a  given  scattering  point  generally  will  be  at  different 
ranges  from  the  source  and  receiver  owing  to  the  horizontal  source-receiver  separation,  and  therefore  a 
different  set  of  rays  will  connect  the  scattering  point  to  the  source  than  to  the  receiver. 


COMPUTER  IMPLEMENTATION 

A  four  stage  sequence  of  computer  programs  was  developed  to  numerically  calculate  boundary 
reverberation.  The  first  three  stages  provide  preliminary  calculations  necessary  for  the  evaluation  of 
Eqs.  (10)  or  (15).  These  equations  are  numerically  evaluated  in  the  fourth  stage.  The  process  is  illus¬ 
trated  by  the  block  diagram  in  Fig.  2. 


Implementation  of  the  four  stage  sequence  of  computer  programs  for  calculating  bistatic  rever¬ 
beration  is  as  follows: 

1.  The  medium  is  described  by  providing  a  sound  speed  profile  and  a  bottom  depth  to  the  com¬ 
puter  program  representing  the  first  stage.  The  sound  speed  profile  consists  of  discrete  depths  for 
which  the  sound  speed  is  specified.  An  output  file  is  produced. 
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STAGES 


RANGE  OF  PROFILES 
SOUND  SPEED  PROFILES 
BOTTOM  DEPTH  vs  RANGE  POINTS 


VELOCITY-FIELD 

CONSTRUCTION 


PROFILE  FILE 


RAY-TRACING 


RAY-TRACE  FILE 


REORGANIZATION 
AND  CONTOUR 
CONSTRUCTION 


CONTOUR  FILES 


REVERBERATION 

CALCULATION 


SOUND  SPEED  AT  DISCRETE  RANGES 
SOUND  SPEED  WITH  RESPECT  TO  DEPTH 
BOTTOM  DEPTH  vs  RANGE  POINTS 


BOTTOM-LOSS  TABLES 
INITIAL  ANGLES 
ITERATION  CONTROLS 


BOUNDARY  ENCOUNTERS  AND  TURNING  POINTS, 

•  INITIAL  ANGLE  •  TRAVEL  TIME 

•  ORDER  •  BOUNDARY  GRAZING  ANGLE 

•  RANGE  •  BOTTOM  LOSS 

•  MAXIMUM  GRAZING  ANGLE 

~|  I  SOURCE  AND  RECEIVER  BEAM  PATTERNS 


ORDER  CONTOURS  (SURFACE  OR  BOTTOM  OR  BOTHI 
TRANSMISSION  LOSS  CALCULATIONS 
STATISTICAL  RAY  AVERAGING  IGPTIONALI 


BOUNDARY  SCATTERING  STRENGTH 

SIGNAL  DURATION 

SOURCE  RECEIVER  SEPARATION 


REVERBERATION  FILES 


PLOTTING  AND 
AVERAGING 


SURFACE  OR  BOTTOM  ENVELOPE 
VERTICAL  DISTRIBUTIONS  OF 
REVERBERATION 


SOURCE  LEVEL 
NOISE  LEVEL 

EFFECTIVE  AZIMUTHAL  BEAMWIDTH 


Fig.  2  —  Four-stage  sequence  of  computer  programs  to  calculate  boundary  reverberation 


2.  A  ray  tracing  program  is  executed  to  calculate  all  boundary  hits  associated  with  rays  that  are 
traced.  The  input  to  this  program  consists  of  the  file  just  created  plus  additional  card  inputs  which 
specify  ray-tracing  parameters  such  as  initial  angles  of  rays  to  be  traced,  various  iteration  controls,  and 
bottom-loss  data.  Bottom-loss  information  is  specified  in  the  form  of  tables  of  grazing  angle  vs  dB  loss. 


The  ray-tracing  program  is  run  twice,  tracing  rays  from  the  source  to  the  boundaries  and  then 
tracing  rays  from  the  receiver  to  the  boundaries.  Each  execution  produces  an  output  data  file  contain¬ 
ing  information  on  the  location  and  character  of  the  ray  reversals  encountered  during  ray  tracing.  By 
character  is  meant  ray  order  (see  below),  initial  angle  of  the  ray,  travel  time,  ray  grazing  angle  at  a 
boundary  encounter,  and  accumulated  bottom  loss  along  the  ray. 

3.  The  ray-tracing  results  are  reorganized,  and  the  transmission  loss  along  each  ray  is  calculated  at 
each  computed  boundary  hit.  Surface  and  bottom-hit  data  are  processed  separately.  Volume  absorption 
and  geometric-spreading  loss  are  combined  with  the  accumulated  bottom  loss  to  determine  the 
transmission  loss  for  each  ray.  At  this  point  the  caustic  correction  is  applied  if  desired.  These 
transmission-loss  values  are  weighted  by  the  source  or  receiver  vertical  beam  pattern.  The  results  of 
the  reorganization  and  the  computed  transmission  losses  are  written  onto  data  files. 
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4.  The  data  files  containing  the  reorganized  ray-tracing  results  are  input  to  a  program  that  calcu¬ 
lates  the  received  bistatic  boundary  reverberation  from  either  the  surface  or  bottom  at  a  discrete  set  of 
times  by  means  of  Eq.  (15).  The  boundary  scattering  strength  is  specified  in  tabular  form  as  a  function 
of  incident  and  backscattered  angles.  Although  not  available  in  the  version  of  the  reverberation  pro¬ 
duction  model  reported  here,  the  authors  have  expanded  the  reverberation  calculation  to  include  esti¬ 
mates  of  both  doppler  shifting  of  surface  reverberation  and  frequency  spreading  of  all  surface  forward 
scattered  returns  in  either  surface  or  bottom  reverberation.  This  extension  of  the  model  will  be 
reported  elsewhere. 

The  basic  reverberation  program  produces  the  calculated  reverberaiior  'velope  vs  time  on 
another  data  file.  In  turn  this  file  is  used  by  a  program  that  plots  the  predicted  .  ve  after  adjusting  for 
source  level  and  noise  level.  The  plotting  program  is  also  capable  of  averagi  several  curves  and  of 
time  averaging.  It  is  also  possible  to  obtain  the  vertical  angular  arrival  structure  trie  reverberation  at  a 
set  of  discrete  times.  Future  inclusion  of  doppler  and  spreading  effects  will  pr  e  the  additional  capa¬ 
bility  of  producing  a  series  of  intensity  vs  frequency  vs  time  plots. 

The  bistatic  calculation  can  be  collapsed  to  compute  "monostatic'  or  "pseudo  monostatic"  rever¬ 
beration  characteristics  by  appropriate  selection  of  input  parameters.  Monostatic  refers  to  the  usual  case 
of  coincident  source  and  receiver.  The  case  where  source  and  receiver  are  separated  in  depth  but  not 
horizontally  is  referred  to  as  pseudo  monostatic.  In  either  of  these  cases,  a  range-dependent  environ¬ 
ment  may  be  specified  by  inputting  to  stage  1  a  sequence  of  sound-speed  profiles  and  bottom  depths  vs 
range.  We  now  describe  each  of  these  four  stages  in  more  detail. 

Sound-Speed  Profile  and  Bottom  Depth  Contour 

Before  the  actual  ray  tracing  can  be  carried  out,  the  sound-speed  field  and  bottom  topography 
must  be  modeled  and  written  onto  a  file.  This  file  serves  as  an  input  to  program  REVRAP.  the  ray¬ 
tracing  program.  The  sound-speed  field  construction  is  accomplished  by  program  PROFIL.  If  necessary 
the  profile  is  extended  linearly  from  the  last  depth  input  to  the  bottom  depth  Program  PLTENV 
creates  Calcomp  plots  of  the  profiles  and  bottom  topography  as  a  function  of  range  from  the  file  gen¬ 
erated  by  PROFIL. 

Mathematical  details  of  the  sound-speed  field  and  bottom  topography  construction  are  in  [4,  p. 
52], 

Ray  Tracing 

Program  REVRAP  reads  the  file  created  by  program  PROFIL  and  traces  rays  using  a  method 
based  on  a  technique  suggested  by  Hudson  Laboratories  [4).  The  procedure  is  an  iterative  one  in  which 
a  ray  is  incremented  from  point  to  point  along  its  path.  This  is  accomplished  by  evaluating  Taylor- 
series  expansions  in  arclength  of  various  ray  parameters  such  as  range,  depth,  travel  time,  and  ray 
angle.  These  expressions  are  derived  from  the  basic  ray  equation: 


d_ 

i  dp 

=  V 

1 

ds 

c(x.  : )  ds 

c(x.  :) 

The  sound  speed  r(.v,  r)  is  assumed  to  be  known  at  every  range  v  and  depth  c  throughout  the  two 
dimensional  medium.  P  is  the  position  vector  to  a  point  on  the  ra>  •  I  s  is  arclength  along  the  rav 

The  ray-tracing  model  has  the  following  features1 

•  Accommodates  multiple  profiles  In'-  possible  horizontal  variations  in  sound  speed 
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•  F prutile  is  defined  by  a  table  of  sound  speed  vs  depth,  with  weighted  parabolic  interpolation 
used  between  specified  depths 

•  Assumes  a  linearly  segmented  ocean  bottom  and  fiat  surface,  with  rays  specularly  reflecting  from 
both  boundaries 

•  Uses  incremental  finite-series  approximations  to  ray  paths 

»  Multiple  bottom-loss  tables  may  be  applied  to  various  range  intervals 

•  Up  to  1000  rays  can  be  traced. 

For  the  purpose  of  estimating  surface  and  bottom  reverberation,  a  set  of  rays  including  the  source 
and  receiver  main  beam  are  traced.  These  rays  are  chosen  so  that  the  resulting  boundary  hits  approxi¬ 
mate  the  true  order  contours  accurately.  If  for  any  reason  it  is  felt  that  an  insufficient  number  of  rays 
have  been  traced  initially,  additional  ray  traces  may  be  performed  and  the  results  consolidated  by  using 
program  MERGE. 

The  rays  are  traced  through  the  medium  which  is  described  by  multiple  sound-speed  profiles,  a 
linearly  segmented  bottom,  and  a  flat  surface.  They  are  traced  one  at  a  time  between  certain  predeter¬ 
mined  ranges,  after  which  accumulated  boundary-hit  and  turning-point  (crest  and  valley)  information  is 
written  onto  a  file.  For  each  ray  reversal,  seven  statistics  are  recorded: 

•  Initial  angle  of  ray 

•  Order  of  the  occurrence 

•  Range  of  the  occurrence 

•  Travel  time 

•  Ray  grazing  angle  at  boundary  encounters  or  the  depth  of  turning  point 

•  Accumulated  bottom  loss 

•  Maximum  grazing  angle  with  which  the  ray  has  encountered  the  surface. 

Specific  details  of  the  ray  tracing  iteration  are  in  [4,  pp.  52,  55-56] 

Reorganization  and  Transmission-Loss  Calculations 

Each  file  created  by  an  execution  of  program  REVRAP  serves  as  an  input  file  to  programs 
PLTCON  and  CONTUR.  Program  PLTCON  provides  a  Calcomp  plot  of  "order  contours"  as  a  function 
of  initial  angle  and  horizontal  range  by  reordering  the  rays  traced  by  program  REVRAP.  (See  explana¬ 
tion  of  order  contours  below.)  Program  CONTUR  also  reorders  the  rays  traced  by  REVRAP  to  form 
order  contours.  In  addition  program  CONTUR  computes  the  full  transmission  loss  for  each  ray  at 
every  boundary  encounter  with  an  optional  caustic  correction  and  creates  an  output  file.  This  output 
file  is  a  necessary  input  file  for  program  BISTRV  which  calculates  reverberation. 

Order  Contours  —  We  first  describe  the  concept  of  order  contours  as  they  relate  to  a  monostatic 
reverberation  calculation.  Consider  the  cyclic  ray  propagated  from  the  point  S-R  in  Fig.  3.  We  adopt 
the  convention  of  numbering  the  ray  crests  or  surface  hits  with  odd  integers  and  using  even  integers  for 
the  valleys  or  bottom  hits.  The  integers  increase  away  from  the  source  and  are  used  to  classify  the  ray 


8 


m jasmwm ~a 


V  _V 


Fig.  3  —  Numbering  of  ray  crests,  ray  valleys,  and  boundary  hits 
illustrated  for  one  ray 


by  the  number  of  oscillations  it  has  made.  Rays  with  the  same  number  of  oscillations  are  said  to  be  of 
the  same  "order."  Order  contours  then  are  derived  from  the  ranges  at  which  rays  encounter  a  boundary. 
Ranges  corresponding  to  ray  reversals  (i.e. ,  ranges  at  which  rays  reverse  their  vertical  direction)  which 
are  not  due  to  boundary  encounters  do  not  define  points  on  order  contours.  Order  contours  are  illus¬ 
trated  in  Fig.  4  where  curves  determined  by  the  boundary  encounters  of  a  given  order  are  plotted  on  an 
initial  source-angle- vs-range  coordinate  system.  A  given  contour  need  not  necessarily  be  a  continuous 
curve  but  may  consist  of  several  disjoint  segments  separated  by  intervals  of  ray  turning  points. 
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Fig.  4  —  Order  contours 

For  the  reverberation  from  a  given  boundary,  only  those  contours  corresponding  to  that  boundary 
are  pertinent;  they  are  odd  for  the  surface,  and  even  for  the  bottom.  For  example,  to  determine  all  ray 
paths  between  the  source-receiver  and  the  bottom  at  range  x,  consider  a  horizontal  slice  at  range  x,  as 
shown  in  Fig.  5.  We  denote  the  three  paths  having  initial  angles  of  by  1,  2,  3,  respectively. 

Then  the  index  set  mix)  of  Eq.  (10)  consists  of  the  integers  1,  2,  and  3.  Therefore,  there  are  32  or  9 
possible  routes  from  the  source-receiver  to  the  bottom  at  range  x  and  back  again. 
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The  order  contours  shown  in  Fig.  6  represent  the  kind  of  bottom  contours  obtained  for  a  typical 
deep  ocean  sound-speed  profile.  The  points  of  maximum  range  on  each  bottom  order  contour 
correspond  to  rays  that  just  graze  the  bottom. 

Figure  7  illustrates  a  smoothed  surface  order  contour  obtained  for  a  typical  deep  ocean  sound- 
speed  profile.  The  two  rays  that  just  graze  the  bottom  (corresponding  to  upward  and  downward  initial 
angles)  at  ranges  XB  and  XB'  determine  points  B  and  B'  on  the  contour.  The  two  rays  that  just  graze 
the  surface  at  ranges  XA  and  XA  occur  at  shallower  initial  angles.  These  rays  contribute  points  A  and  A ' 
on  the  contour.  Thus  rays  corresponding  to  initial  angles  in  the  intervals  (9#,  6A •)  and  (0A,  oB )  hit 
the  surface  but  not  the  bottom.  Rays  providing  the  significant  contributions  to  surface  reverberation  at 
long  ranges  have  initial  angles  in  the  intervals  defined  above.  Rays  having  steeper  initial  angles  hit  the 
bottom,  thereby  suffering  bottom  loss.  The  cumulative  loss  associated  with  several  bottom  hits  causes 
these  rays  to  have  less  effect  on  surface  reverberation  at  long  ranges  where  the  higher  order  contours 
are  relevant. 
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Fig.  7  —  Smoothed  surface  order  contour  illustrating  caustic  behavior 


The  bistatic  and  pseudo  monostatic  cases  require  two  ray  traces,  one  for  the  source  and  one  for 
the  receiver  Lach  boundary  point  now  corresponds  to  two  distinct  distances,  source-to-boundary  point 
distance  and  boundary  point-to-receiver  distance.  Thus,  for  each  boundary  point,  the  values  assumed 
by  index  i  of  Eq  (15)  are  determined  by  the  properties  of  the  source  ray  trace  while  the  values 
assumed  by  index  j  of  Eq.  (15)  depend  upon  the  receiver  ray  trace.  Two  sets  of  order  contours  are 
required  for  each  boundary  (surface  and  bottom). 

Transmission-Loss  Calculations  —  Initially,  an  estimate  of  the  ray  intensity  is  made  at  each  location 
on  the  contour  determined  by  the  ray  tracing.  Transmission  loss  consists  of  bottom  loss,  volume 
absorption,  and  geometric-spreading  loss.  Bottom  loss  is  computed  during  ray  tracing;  whereas,  the 
remaining  losses  are  calculated  after  the  order  contours  are  constructed.  Volume  absorption  is 
estimated  by  using  Thorp's  equation  [4],  Geometric-spreading  loss  under  our  assumption  of  azimuthal 
symmetry  is  given  by  [5]: 
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Ux I  -  4^|i 
a*  cos  H  dtf 


where  a  —  unit  reference  distance,  x  =  horizontal  range,  *  initial  source  angle  of  ray,  and  y=angle 
the  ray  makes  with  the  horizontal  at  range  x. 


The  derivative  dx/d#  in  Eq.  (17)  is  approximated  numerically  at  each  point  (x,,  W,)  on  a  contour. 
If  0,_|,  0„  0,  +  i  are  successive  initial  angles  of  rays  that  are  traced  on  this  contour  and  x,_|,  x,,  x,  +  |  are 
the  corresponding  ranges,  then 


Ax,  A  R,  A 


A  R,-t  A  R, 


A#,  A R  t),  -  A R  «)  +  |  -  0,  ' 


where 

A/?,_i  —  |x(  -  x,_j |,  (19) 

A /?,  -  |x/+i  -  x,l,  and 
A  R  -  A/?,_i  +  A/?^ , 

are  used  to  estimate  |9x/90|.  Transmission  losses  between  computed  boundary  encounters  will  be 
found  by  interpolating  linearly  in  range  between  the  estimates  at  the  computed  boundary  hits. 

In  using  a  ray  approach  to  the  problem  of  sound  propagation,  a  serious  deficiency  arises  when 
attempting  to  calculate  geometric-spreading  loss  in  near-caustic  regions,  that  is,  near  locations  where 
adjacent  ray  paths  cross.  In  this  region  geometric  spreading  tends  to  zero,  erroneously  implying  infinite 
intensity  at  the  point  of  intersection.  This  phenomenon  occurs  when  the  derivative  6x/9«  is  equal  to 
zero,  or,  as  related  to  order  contours,  where  the  contour  of  order  n,  x  =  /„(0),  has  a  stationary  value. 
The  surface  order  contour  shown  in  Fig.  7  has  caustics  at  points  C  and  C'. 

Also,  as  rays  are  traced  to  longer  ranges  from  the  source,  the  corresponding  order  contours  tend 
to  become  irregular,  owing  to  accumulated  computation  errors  and  numerical  approximations.  As  a 
result,  computed  transmission  loss  for  the  higher-order  contours  becomes  increasingly  erratic  and 
uncertain. 

To  handle  these  problems,  a  statistical  ray-averaging  procedure  has  been  developed  that  operates 
on  the  order  contour  and  effectively  smooths  irregularities  in  the  propagation-loss  predictions.  It 
involves  unfolding  the  order  contour  and  providing  an  avenue  around  zeros  in  dx/do.  The  procedure 
leads  to  an  expression  that  is  similar  to  the  result  obtained  by  Hardy  [51.  Whereas  Hardy’s  results  call 
for  ray  intensity  averaging  as  a  function  of  depth,  the  method  used  in  program  CONTUR,  developed  by 
Palmer  [61,  performs  ray  intensity  averaging  as  a  function  of  range. 

The  presence  of  caustics  corresponding  to  rays  having  initial  angles  within  the  intervals  defined  by 
the  surface  grazing  ray  angles  and  the  bottom  grazing  ray  angles  implies  that,  for  surface  reverberation, 
the  greatest  variations  in  ray  intensity  occur  in  those  intervals.  Thus  the  initial  angles  of  rays  that  con¬ 
tribute  to  the  corresponding  significant  portions  of  the  order  contours  must  be  sampled  densely  enough 
during  ray-tracing  to  ensure  accurate  intensity  calculations. 

The  total  transmission  loss  at  the  surface  or  bottom  boundary  as  a  function  of  horizontal  range 
from  the  source  or  receiver  can  be  computed  by  using  program  TOTLOS.  The  program  produces  a 
Calcomp  plot  of  transmission  loss  vs  range. 
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Reverberation  Calculation 

Acoustic  reverberation  from  either  the  ocean  surface  or  bottom  is  calculated  in  program  BISTRV. 
The  calculation  is  performed  by  numerically  evaluating  Eq.  (15).  Program  BISTRV  requires  two  data 
files  as  input,  each  created  by  an  execution  of  program  CONTUR.  The  two  files  correspond  to  ray 
information  from  the  source  and  from  the  receiver  with  respect  to  either  the  surface  or  bottom;  i.e., 
one  execution  of  program  BISTRV  calculates  reverberation  from  either  the  ocean  surface  or  ocean  bot¬ 
tom. 


The  choice  of  a  coordinate  system  for  the  evaluation  of  Eq.  (15)  depends  in  part  on  the  underly¬ 
ing  ray-trace  programs.  Since  the  ray-trace  information  such  as  transmission  loss  and  travel  time  is 
ordered  by  range  from  the  source  or  receiver,  it  is  advantageous  to  use  a  coordinate  system  ( p ,  r )  in 
which  p  is  the  horizontal  range  along  the  propagation  path  out  from  the  source,  and  r  is  the  horizontal 
range  along  the  return  path  back  to  the  receiver.  The  ordered  pair  ( p ,  r)  constitutes  a  biradial  coordi¬ 
nate  system.  The  transformation  of  Eq.  (15)  to  the  biradial  coordinate  system  yields  (see  the  appen¬ 
dix): 


x  a(9!,  <f>„  9j,  <t>j)  [Sh2(p2  +  r2)_P{p2*Pr2)2_  16A4]l/2 


(20) 


where  2 h  is  the  source-receiver  horizontal  separation.  The  source-receiver  geometry  and  the  region  of 
integration  are  shown  in  Fig.  8. 


v 


Fig.  8  —  Source-receiver  boundary  plane  geometry  for  bistatic 
situation  in  terms  of  both  Cartesian  and  biradial  coordinate 
systems 


Equation  (20)  is  evaluated  numerically  by  incrementing  the  ranges  p  and  r  in  equal  step  sizes 
from  the  source  and  receiver  respectively.  At  present,  the  program  assumes  that  the  horizontal  com¬ 
ponent  of  the  beam  patterns  of  both  source  and  receiver  are  omnidirectional.  This  decreases  the  com¬ 
putation  time  significantly.  Equation  (20)  now  becomes 
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The  value  R(t„ )  is  defined  to  be  the  reverberation  density  per  A t  occurring  in  the  time  interval 
[t„,  rn+il,  where  A t  =  rn+!  -  r„,  is  the  time  increment.  The  values  t0,  iN  and  A t  are  input  parameters 
of  program  BISTRV.  Values  of  the  various  ray  parameters  at  each  boundary  point  (pk,  /•,)  are  obtained 
by  linearly  interpolating  between  points  on  the  appropriate  order  contours. 

The  scattering  strength  a  is  approximated  by  using  a  backscattering  table  that  is  input  by  the  user. 
This  table  defines  backscattering  as  a  function  of  grazing  angle.  The  scattering  value  used  for  <t  is  that 
value  in  the  table  that  corresponds  to  the  average  of  the  grazing  angles  and  of  the  rays  along  paths 
/  and  j  at  the  boundary  point  ( pk ,  r,).  The  scattering  strength  it  is  assumed  to  be  independent  of  the 
azimuthal  angles  <£,  and 

Each  term  of  the  double  summation  in  braces  (/,  j  summation  in  Eq.  (21))  represents  the  contri¬ 
bution  to  reverberation  due  to  a  particular  reverberant  path  i,  j  from  source  to  boundary  area  (defined 
by  ( pk ,  r t))  to  receiver.  As  each  contribution  is  calculated,  the  time  interval  over  which  it  will  be 
received  is  computed.  The  contribution  is  then  added  to  those  cells  of  the  time  array 
R(t„),  n  =  0, 1 . )V,  for  those  time  values  t„  falling  within  the  computed  time  interval  of  recep¬ 

tion.  Thus  each  calculated  contribution  is  added  to  one  or  more  of  the  cells  of  array  R(tn)  (  provided 
of  course  that  the  time  interval  of  reception  does  not  lie  outside  of  the  time  interval  of  interest 
\t0,  fjv  +  A/1).  If  the  time  interval  of  reception  of  a  contribution  only  partially  spans  a  time  interval 
[/„,  t„  +  Ar]  of  the  R(tn)  array,  then  the  contribution  is  reduced  by  that  fraction  of  A t  before  being 
added  to  that  cell.  After  all  the  reverberation  contributions  have  been  calculated,  each  value  of  the 
R  (/„)  array  is  divided  by  A  t  to  obtain  reverberation  density  per  unit  time. 

The  above  calculation  is  valid  for  the  true  bistatic  case  where  the  source  and  receiver  are 
separated  in  range,  only  if  the  environment  is  range  independent;  i.e.,  the  ocean  bottom  is  flat  and  one 
sound-speed  profile  is  defined.  This  is  because  of  the  use  of  two-dimensional  (range- vs-depth)  ray  trac¬ 
ing  and  assumed  azimuthal  symmetry  of  acoustic  propagation  from  both  the  source  and  receiver.  When 
the  source  and  receiver  are  not  coincident  in  range,  concurrent  azimuthal  symmetry  from  the  source 
and  receiver  is  possible  only  in  a  range  independent  environment.  However,  when  the  source  and 
receiver  coincide  in  range,  the  environment  may  be  range  dependent. 

Program  BISTRV  is  also  capable  of  providing  the  vertical  distribution  of  reverberation  at  specific 
vertical  angles  and  times  input  by  the  user.  The  angles  at  which  reverberation  values  are  provided  are 
with  respect  to  the  receiver.  The  reverberation  contribution  at  angle  0  at  time  t  is  that  portion  of  the 
total  reverberation  density  corresponding  to  rays  having  acceptance  angles  in  the  angle  interval  of  size 
AO  (input)  centered  at  0,  divided  by  AO.  This  yields  reverberation  density  per  unit  time  per  unit  angle. 

A  condensed  form  of  Eq.  (21)  is  evaluated  for  pseudo  monostatic  and  monostatic  reverberation 
(see  description  of  program  BISTRV).  The  calculation  assumes  full  360°  azimuthal  omnidirectionality. 
However,  the  user  can  provide  for  horizontal  beam  widths  by  adjusting  the  value  of  source  level  input 
in  program  AVREVB  which  is  discussed  in  the  next  paragraph. 

The  reverberation  values  are  output  on  a  file  which  is  used  by  program  AVREVB.  This  program 
plots  a  curve  of  reverberation  vs  time.  Values  of  source  level  (possibly  adjusted  for  azimuthal  beam 
width  for  a  monostatic  calculation)  and  noise  level  are  input  parameters  of  program  AVREVB.  The 
program  is  also  capable  of  averaging  several  reverberation  curves  and  of  time  averaging. 


DESCRIPTION  OF  PROGRAMS 

The  computer  programs  described  in  this  report  are  currently  operational  on  the  Texas  Instru¬ 
ments  Advanced  Scientific  Computer  (ASC)  at  the  Naval  Research  Laboratory.  This  section  describes 
the  purpose  and  usage  of  each  program  involved  in  the  sequence  of  calculations  and  plottings. 
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Program  Name:  PROFIL 

Purpose:  Calculate  velocity  gradients  and  curvatures  as  functions  of  depth  from  input  velocity  profiles. 

Create  an  output  file  containing  the  velocity  profiles,  their  calculated  derivatives,  and  a 


description  of  the  bottom  topography. 

Usage: 

Input 

1.  Punched  Cards 

Format 

Card  #1  TITLE 

TITLE  —  80  column  title  for  printed  output. 

(20A4) 

Card  #2  NECUR 

(15) 

NECUR—  flag  for  earth’s  curvature  correction. 

—  1  if  curvature  correction  is  desired, 

—  0  otherwise. 


Card  #3  ND1,  NMO,  NDY,  NYR,  NPRF 
ND1  —  arbitrary  identification  number. 

NMO  —  month 

NDY  —  day  arbitrary  date. 

NYR  —  year 

NPRF  —  number  of  velocity  profiles  to  be  input  (NPRF>2). 

Card  #4  NBPTS,  BOTMIN,  BOTMAX,  RMAX 

NBPTS  —  number  of  bottom  points  used  to  describe  a  linearly 
segmented  bottom  (2  <  NBPTS  <  500). 

BOTMIN  —  miminum  bottom  depth  (m). 

BOTMAX  —  maximum  bottom  depth  (m). 

RMAX  —  maximum  range  (nmi)  of  bottom  points. 


Type  #5  (BTRNGE(I),  BTDPTH(I),  I  -  1,  NBPTS)  (2F10.4) 

BTRNGE  —  range  (nmi)  of  bottom  point. 

BTDPTH  —  depth  (m)  of  bottom. 

Note:  NBPTS  cards  of  type  #5  are  read. 

Type  #6  R  (F10.1) 

R  —  horizontal  range  (nmi)  at  which  the  velocity  profile 
to  be  read  in  is  defined. 


Type  #7  (D(I),  V(I),  I  -  1,M)  (10F8.1) 

D  —  depth  (m). 

V  —  velocity  (m/s). 

Note:  Type  #7  is  repeated  until  M  nonzero  velocities  are  read  in 
at  a  rate  of  up  to  5  per  card.  (M<500!> 


(515) 


(I5,F15.4,2F10.4) 
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1.  Punched  Cards  Format 

Type  #8  —  1  in  columns  14  and  15  (115) 

signals  that  there  are  no  more  points  in  this 
velocity  profile. 

Note:  Types  #6-8  are  repeated  until  NPRF  profiles  are  read  in. 

For  a  single  profile  environment  the  profile  is  read  in  twice, 
once  at  zero  range  and  once  at  maximum  range. 

Type  #9  —  7  in  columns  14  and  15  (115) 

signals  that  there  are  no  more  profiles  to  be  input. 

Subroutines  Called:  DERIV,  PARAM 


Output 


1.  Printed:  Velocity  profiles  with  gradients  and  curvatures. 

2.  File:  FT01F001.  Catalogue  this  file  for  use  as  input  for  programs  REVRAP  and  PLTENV. 


Program  Name:  PLTENV 

Purpose:  Produce  a  Calcomp  plot  of  velocity  profiles  and  bottom  topography  from  a  file  created  by  pro¬ 
gram  PROFIL. 


Usage: 


Input 

1.  Punched  Cards  Format 

Card  #1  RM,  DM,  XAXIS,  YAXIS,  VTR  (5F10.3) 

RM  —  maximum  range  (nmi)  of  plot. 

DM  —  maximum  depth  (m)  of  plot. 

XAXIS  —  length  of  x-axis  of  plot  in  inches. 

YAXIS  —  length  of  y-axis  of  plot  in  inches. 

VTR  —  velocity  to  range  scale  for  plotting.  VTR  equals  the 

length  in  inches  that  represents  10  m/s  for  the  velocity 
profile  curves. 

2.  Files:  Assign  output  file  from  program  PROFIL  to  FT01F001. 

Subroutines  Called:  LABEL,  PLOT  AX,  PLOTPT 
Output 

Plot:  A  YAXIS  by  XAXIS  Calcomp  plot  is  generated.  The  plot  includes  bottom  depth  (m) 
versus  range  (nmi)  and  the  velocity  (m/s)  vs  depth  (m)  profiles  located  at  their  corresponding 
ranges  (nmi). 
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Program  Name:  REVRAP 

Purpose:  Trace  rays  using  iterative  technique.  Calculate  accumulative  bottom  loss  and  part  of  the 
geometric  spreading  loss  formula  at  each  boundary  hit  or  refractive  reversal.  Calculate 
travel  time  and  horizontal  range  along  with  either  the  grazing  angle  for  a  boundary  hit  or  the 
reversal  depth  for  a  ray  reversal  caused  solely  by  refraction.  Obtain  the  maximum  surface 
grazing  angle  for  each  ray.  Generate  output  file  containing  ray  information  at  reversal  points 
(refractions  and  boundary  hits). 

Usage 


1 .  Punched  Cards 


Format 


Card  #1  TITCAR 

TITCAR  —  80  column  title  for  printed  output. 


(20A4) 


Card  #2  RNGINI,  DPINIT,  DELMAX,  DELMIN,  VFAEPS, 
DELSMX,  SURFDE,  TIMEMA 


(8F10.3) 


RNGINI 

DPINIT 

DELMAX 

DELMIN 

VFAEPS 

DELSMX 


SURFDE 


TIMEMA 


range  of  source  (nmi)  (usually  zero), 
depth  of  source  (m). 

first  trial  step  size  (m)  used  in  predicting  a  new 
ray  path  position  (e.g.,  1000  m). 
minimum  allowable  step  size  (e.g.,  20  m). 
velocity  field  accuracy  (e.g.,  0.2  m/s), 
maximum  allowable  change  in  the  sine  of  the  ray 
angle  (angle  from  horizontal)  between  predicted 
positions  (e.g.,  0.02). 

vertical  distance  (m)  from  surface  or  bottom  within 
which  a  ray  is  to  be  considered  as  hitting  that 
boundary  (e.g.,  0.05  m). 

travel  time  (s)  after  which  a  ray  will  be  terminated. 


Card  #3  TLOSSM 


(F10.3) 


TLOSSM  —  intermediate  transmission  loss  value  (dB)  after  which 
a  ray  will  be  terminated  (e.g.,  200  dB). 

Card  #4  NGRPS,  NRGRPS,  KTAPE,  MAXORD 

NGRPS  —  number  of  cards  of  type  #5  to  be  read  in  to  define 

all  the  initial  angles  of  the  rays  to  be  traced. 

NRGRPS  —  rectification  range  identifier. 

if  positive  —  number  of  cards  of  type  #6  to  be  read 
in  to  define  "rectification"  ranges. 


if  negative  —  negative  of  number  of  rectification 
ranges  to  be  read  in  on  cards  of  type  #7  at  a  rate 
of  8  ranges  per  card. 
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1.  Punched  Cards  Format 

Note:  As  each  rectification  range  is  reached,  ray  information  obtained  up  to  that  range  is 
printed  and  written  on  an  output  file.  Any  range  at  which  a  velocity  profile  has  been 
defined  by  program  PROFIL  is  considered  an  additional  rectification  range  to  those 
defined  by  card  types  #6  or  #7. 

’  KTAPE  -  file  identifier 

if  positive  —  logical  unit  of  output  file  if  desired 
(KTAPE  ^  2  but  *  5,  6,  or  7). 
if  nonpositive  —  no  output  file  is  generated. 

Note:  Output  file  will  be  assigned  to  FTONFOOl  where  N  =  KTAPE. 


MAXORD  —  maximum  allowable  boundary  encounter  order  to  which 
each  ray  is  traced. 

Type  #5  ANGINC,  ANGIN1,  ANGFIN  (3F10.3) 

ANGINC  —  angle  increment  (degrees)  of  group. 

ANGINI  —  initial  angle  (degrees)  in  group. 

ANGFIN  —  final  angle  (degrees)  in  group. 

Note:  The  positive  Z  axis  is  directed  downward;  i.e.,  positive 

angles  describe  rays  whose  initial  direction  is  towards 
the  bottom.  NGRPS  cards  of  type  #5  are  input  with  two 
requirements: 

1.  Angles  must  be  defined  in  increasing  order 
from  negative  angles  to  positive  angles, 


2.  An  angle  can  be  defined  in  at  most  one  group, 


Type  #6 

RGINC 

RGI 

RGF 

Note: 


Type  #7 


i.e.,  no  overlapping  by  angle  groups  is  permitted. 

The  maximum  number  of  angles  allowed  is  1000. 

RGINC,  RGI,  RGF  (3F10.3) 

—  rectification  range  increment  (nmi). 

—  initial  rectification  range  (nmi). 

—  final  rectification  range  (nmi). 

NRGRPS  cards  of  type  #6  are  input  only  if  NRGRPS  is  positive.  The  first 
rectification  range  defined  should  be  zero.  A  rectification  range  can  be  defined  by  at 
most  one  card  of  type  #6,  i.e.,  no  overlapping.  A  maximum  of  8000  total  reversals 
are  allowed  in  each  rectification  interval. 

(STRRNG(I),  I  -  1,  NRGRPS)  (8F10.3) 


STRRNG(I)  -  Ith  rectification  range  (nmi).  (NRGRPS  <  300) 


NRL  REPORT  8721 


1.  Punched  Cards  Format 

Note:  (-NRGRPS)  ranges  are  input  at  a  rate  of  8  per  card  only  if  NRGRPS  is  negative.  The 
first  rectification  range  defined  should  be  zero. 


Card  #8 
NINT 

NTAB 


NINT,  NTAB,  (NINTAB(I),  1=1,  NTAB) 

—  number  of  range  intervals  for  which  different 
bottom  loss  tables  are  desired  (1)  <  NINT<10). 

—  number  of  distinct  bottom  loss  tables  (1<NTAB<4). 


(615) 


NINTAB(I)  —  number  of  loss  values  in  Ith  table  (NINTAB(I<  100). 


Type  #9  GRANG(I.J),  DBLOSS(I,J)  (2F10.3) 

GRANG—  grazing  angle  (degrees)  (0  <  GRANG  <  90). 

DBLOSS—  corresponding  bottom  loss  in  positive  dB. 

Note:  First  NINTAB(l)  cards  of  type  #9  are  input  to  define  loss  table  #1,  then  NINTAB(2) 
cards  of  type  #9  are  input  to  define  table  #2,  etc.,  until  NTAB  bottom  loss  tables  are 
input. 

Type  #10  KTAB(l),  ROFTAB(l)  (I2,F8.1) 

KTAB(l)  —  bottom  loss  table  number  to  be  used  in  range  interval 

I  (  1  <  1  <  NINT).  Table  number  is  determined  by  the 
order  in  which  tables  are  read  in. 

ROFTAB(f)  —  maximum  range  (nmi)  of  1th  range  interval. 


Note:  NINT  cards  of  type  #10  are  read,  i.e.,  one  for  each  range  interval  the  Ith  range  inter¬ 
val  is  given  by  /,  —  (R^^Rj),  1  <  /  <  NINT,  where  R0  =  0  and  /?,  = 
ROFTABU+l). 

2.  Files:  Assign  output  file  from  program  PROFIL  to  FT01F001. 


Subroutines  Called  -  HITBOT,  HITSUR,  ITRAT,  RECORD 


Output 
1.  Printed: 

a.  Listing  of  input  parameters. 

b.  A  table  containing  ray  information  at  all  reversal  points  (crests,  valleys, 
surface  hits,  and  bottom  hits).  The  information  is  printed  out  for  increasing 
initial  angle  (angle  at  source)  for  each  retification  interval.  A  description 

of  the  nine  columns  of  the  table  is  as  follows: 
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Label 

REVERSAL 

ANGLE 

CODE 


RANGE  (nmi) 


Description 

A  running  total  of  the  number  of  reversals 
encountered  by  the  program. 

Ray  angle  in  degrees  at  the  source 

An  eight-digit  hexidecimal  word  interpreted  as 
three  distinct  numbers,  AAABBBBC  where 
AAA  is  the  order  of  the  reversal  of  the 
ray  being  described,  BBBB  is  the  ray  number 
which  is  in  one-to-one  correspondence  with 
initial  angle,  C  is  the  type  of  ray  reversal 
where 

1  —  surface  hit 

2  —  bottom  hit 

3  —  crest 

4  —  valley 

Horizontal  range  of  the  reversal  point 
in  nautical  miles. 


GRZ  ANG  OR  DEPTH  Grazing  angle  of  the  ray  in  degrees  at 

the  surface  or  bottom  or  the  depth  in 
meters  of  the  ray  at  its  reversal  point. 

TIME  Travel  time  in  seconds  for  the  ray  to 

reach  the  reversal  point. 

TR  LOSS  Intermediate  transmission  loss  value 

(intensity  units)  obtained  by  evaluating 
the  transmission  loss  formula  (Eq. ( 17))  without 
including  the  dx/d9  factor.  This  factor 
is  calculated  in  program  CONTUR. 

DB  Same  transmission  loss  value  in  dB  units. 

MAX  SURFACE  GRZ  The  maximum  grazing  angle  with  which 

the  ray  has  encountered  the  surface. 

A  value  of  1000  indicates  that  the 
ray  has  not  yet  reflected  from  the 
surface. 


File:  FTONFOOl  where  N  —  KTAPE  (Card  #4).  Catalogue  this  file  for  use  as  input  file  for  pro¬ 

grams  PLTCON,  CONTUR,  and  possibly  program  MERGE. 

Note:  Rays  are  traced  one  at  a  time  from  the  source  to  the  first  rectification  range  (input 
ranges  plus  ranges  of  sound  speed  profiles).  Next  reversal  information  is  printed  and 
written  onto  the  output  file.  The  process  is  repeated  to  each  rectifica,  n  range 
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Program  Name:  MERGE 


Purpose:  Merge  two  or  more  ray  tracing  files  created  by  program  REVRAP  into  a  new  single  ray-tracing 
file  which  can  then  be  used  c«  input  for  programs  PLTCON  and  CONTUR.  Program 
MERGE  enables  the  user  to  increase  the  number  of  rays  traced  from  those  originally  traced 
by  REVRAP  by  executing  REVRAP  again.  The  additional  executions  of  REVRAP  are  used 
to  trace  rays  with  initial  angles  not  traced  by  previous  REVRAP  executions.  This  eliminates 
the  need  to  trace  any  ray  more  than  once. 


Format 


Card  #1  KTAPE,  NOTAPE,  NTAPE  (315) 

KTAPE  —  logical  unit  number  of  main  file. 

NOTAPE  —  number  of  input  files  to  be  merged  (each  file  is  the 
output  file  of  a  REVRAP  execution). 

NTAPE  —  logical  unit  number  of  merged  output  file. 


Type  #2  MTAPE,  NUMIN  (215) 

MTAPE  —  logical  unit  number  of  a  file  to  be  merged  (MTAPE  ^ 5,6,7). 

NUMIN  —  the  number  of  rays  from  MTAPE  file  to  be  merged  with  those 

rays  present  on  file  KTAPE. 

Type  #3  KAY(I),  1  =  1,  NUMIN  (1615) 

KAY  —  Ray  numbers  of  rays  to  be  merged  with  rays  present  on  the 

main  file  (KTAPE).  These  numbers  are  identified  by 
examining  the  hexidecimal  words  labeled  CODE  in  the 
printed  output  of  that  REVRAP  execution  that  created 
the  file  now  assigned  to  logical  unit  number  MTAPE.  Use 
the  decimal  equivalent  of  the  ray  number  as  defined  in 
CODE  (i.e.,  do  not  input  the  ray  number  as  a  hexidecimal 
number).  Ray  numbers  must  be  in  increasing  order,  i.e., 

KAYO+  1)  >  KAY(I). 

Note:  For  each  file  to  be  merged,  one  card  of  type  #2  followed  by  one  or  more  cards  of 
type  #3  are  required.  It  is  not  necessary  to  merge  every  ray  traced  by  an  auxiliary 
REVRAP  execution.  Instead  type  #3  cards  enable  the  user  to  merg  only  those  rays 
that  are  wanted.  Of  course,  the  new  file,  NTAPE,  will  contain  all  rays  that  are 
present  on  the  original  file  KTAPE. 

2.  Files:  Assign  output  file  from  program  REVRAP  which  is  interpreted  as  main  file  to 
FT0KF001 ,  where  K  =  KTAPE. 

Assign  auxiliary  output  files  created  by  additional  executions  of  program  REVRAP  to 
FT0MF001,  where  M  =  MTAPE.  Use  a  different  value  of  MTAPE  for  each  auxiliary 
file.  (M  *  5,  6,  7). 
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Subroutines  Called:  There  are  no  subroutines  called. 


1.  Printed:  Information  describing  the  execution  of  program  MERGE  is  provided  including  the 

angles  for  which  ray  information  exists  on  the  new  (output)  file. 

2.  File:  FTONFOOl,  where  N  =  NTAPE.  Catalogue  this  file  for  use  as  input  file  for 

programs  PLTCON  and  CONTUR. 


Program  Name:  PLTCON 

Purpose:  Produce  Calcomp  plots  of  surface  and/or  bottom  order  contours  from  the  output  file  generat¬ 
ed  by  program  REVRAP. 


1.  Punched  Cards  Format 

Card  #1  (KTITLE(I),  1  =  3,  20)  (18A4) 

KTITLE  —  72  column  title  for  plot.  (The  program  defines  the  first  word  of  the  title  as 

either  SURFACE  or  BOTTOM). 

Card  #2  KTAPE,  (MINORD(I),  MAXORD(I),  1=1,  2)  (515) 

KTAPE  —  logical  unit  number  of  input  file. 

MINORD(l)  —  minimum  order  of  surface  contour  to  be  plotted. 

MAXORD(l)  —  maximum  order  of  surface  contour  to  be 
plotted  (can  be  0). 

MINORD(2)  —  minimum  order  of  bottom  contour  to  be  plotted. 

MAXORD(2)  —  maximum  order  of  bottom  contour  to  be  plotted 
(can  be  0). 

Card  #3  AG1,  AG2,  YAXIS,  RG1 ,  RG2,  XAXIS  (6F10.2) 

AG1  —  minimum  y-coordinate  (initial  angle  —  degree). 

AG2  —  maximum  y  coordinate  (final  angle  —  degree). 

YAXIS  —  length  of  y-axis  in  inches  (  ^  10  in.). 

RG1  —  minimum  x-coordinate  (minimum  range  — nmi). 

RG2  —  maximum  x-coordinate  (maximum  range— nmi). 

XAXIS  —  length  of  x-axis  in  inches. 

2.  Files:  Assign  output  file  from  program  REVRAP  to  file  FTONFOOl  where  N  =  KTAPE  (card 

#2). 


Note.  Present  dimensioning  allows  a  maximum  of  8000  ray  reversals  to  be  input. 


mm 

m 
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Subroutines  Called:  ORDPLT.  PLOTAX.  PLOTPT 

.  *« 

b 

Output 

• 

f\^ 

1.  Printed: 

- 

+m'  -‘ 

Order  contours 

consisting  of  rays  of  the  same  order  encountering  a  boundary  are  listed,  first  all 

*  *  *„* 

surface  contours  then  all  bottom  contours.  Each  contour  is  comprised  of  one  or  more  mono- 

* .  - 

as 

tonic  segments. 

i.e.,  segments  where  the  horizontal  range  of  rays  encountering  the  boundary 

either  increases 

or  decreases  as  a  function  of  initial  angle.  For  each  order  contour,  rays  are 

FT; 

grouped  by  these  monotonic  segments  and  listed  in  the  order  of  increasing  horizontal  range. 

.  < 

\'«\ 

V- 

r>. 

The  segments  are  ordered  by  increasing  initial  angle. 

■_v 

f.’~\ 

■••:••• 

2.  Plots: 

Calcomp  plots  of  order  contours  corresponding  to  surface  and/or  bottom  boundary  encounters 

•1 

’7T 

are  generated. 

The  contours  are  plotted  as  functions  of  initial  angle  (y-axis)  versus  horizontal 

range  (x-axis). 

Surface  contours  and  bottom  contours  are  plotted  on  separate  sets  of  axes. 

s- 

* 

m 

Program  Name:  CONTUR 

v 

55 

Vi. 

r 

Purpose:  Read  output  file  created  by  program  REVRAP,  create  order  contours  and  complete  transmis- 

sion  loss  calculations  by  computing  two  types  of  transmission  loss.  Create  an  output  file 

fcw 

containing  ray  information  including  one  type  of  transmission  loss  calculation  (corrected  for 

.  * 

«E 

source  or  receiver  beam  pattern).  In  contrast  to  program  PLTCON,  each  execution  of  pro- 

j 

fflp 

gram  CONTUR  processes  either  surface  or  bottom  boundary  hits  but  not  both. 

V 

Usage 

to 

If 

Input 

kc- 

ss 

1.  Punched  Cards 

Format 

961 

aw 

Card  #1  TITLE 

(20A4) 

»: 

y/i 

TITLE  —  title  card  for  printed  output. 

■-> 

(515) 

;>«; 

Card  #2  KORD, 

KFRQ,  KTAPE,  LTAPE,  ITL 

v-; 

KORD 

—  type  of  order  contours  to  be  processed 

1 

Jk  ■  .* 

=  1  for  surface  contours 

pi 

?? 

=  2  for  bottom  contours 

['t‘* 

KFRQ 

—  frequency  of  source  in  Hz.  (Used  in  Thorpe’s  equation.) 

KTAPE 

—  logical  unit  number  of  input  file. 

■> 

,  *.  , 

LTAPE 

—  logical  unit  number  of  output  file. 

2 

as 

ITL 

—  flag  for  type  of  transmission  loss  calculation  to  be 
output  on  file. 

• 

*.  .% 

=  1  minimum  of:  a.  statistical  averaging  of  ray 

'■  /» 

.%*% 

intensities,  b.  ray  bundle  estimates  (i.e.,  value 

V  *. 

whi.-h  predicts  the  larger  transmission  loss  is  chosen). 

> 

E 

« 

=  2  ray  bundle  estimates. 

•j 

& 
•  ■*  « 

'%•*<< 

23 

V 

s 

•  v* 
-* F » 

* 

;:c 

*  *  1  .  "  ,  '  ^  .  t  .  •  .  .  .  ■  ^  1  _  * .  '  1  ^  ^  ^  1  *  *  «  *  §  %  |  ^  • 
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1.  Punched  Cards 

Type  #3  (BP(I),  1  =  1,181) 


Format 

(8F10.0) 


BP  —  the  vertical  beam  pattern  of  the  source  or  receiver.  A 

total  of  181  values  are  read  in  by  means  of  23  input  cards.  The 
pattern  is  to  be  specified  from  -90°  (towards  surface)  to 
90°  (towards  bottom)  in  1°  steps.  Values  are  in  dB  units 
with  all  values  less  than  or  equal  to  zero. 

Note:  If  an  Ominidirectional  pattern  with  a  value  of  zero  dB  at  all  angles  is  desired,  then  1 
card  of  type  #3  with  BP(1)  =  99.  may  be  read  in  instead  of  23  blank  cards. 

2.  Files:  Assign  output  file  from  program  REVRAP  to  file  FTONFOOl  where  N  =  KTAPE. 

Note:  Ray  reversals  are  input  in  blocks  corresponding  to  rectification  intervals.  (See  note  fol¬ 
lowing  the  description  of  Card  #4  in  program  REVRAP.)  All  ray  reversals  except 
those  corresponing  to  type  KORD  (Card  #2)  are  deleted  from  each  block  before  the 
next  block  is  input.  The  total  number  of  ray  reversals  present  in  CONTUR  is  at  one 
point  as  large  as  the  number  of  ray  reversals  of  type  KORD  determined  previously  plus 
the  number  of  ray  reversals  in  the  last  rectification  range.  This  total  may  not  exceed 
8000  with  present  dimensioning. 

Subroutines  Called:  FORMOC 

Output 

1.  Printed: 

a.  Information  based  on  input  parameters 

b.  Order  contours  consisting  of  rays  of  the  same  order  encountering  the 
requested  boundary  are  listed.  The  listings  are  similar  to  those  found 
in  the  printed  output  of  program  PLTCON. 

2.  Files:  FTONFOOl  where  N  =  LTAPE  (Card  #2).  Catalogue  this  file  for  use  as  input  file  for 

programs  TOTLOS  and  BISTRV. 

Program  Name:  TOTLOS 

Purpose:  Compute  transmission  loss  (corrected  for  source  or  receiver  beam  pattern)  at  a  boundary 
(surface  or  bottom)  as  a  function  of  range.  Produce  a  Calcomp  plot  of  transmission  loss  vs 
range. 

Usage: 


1.  Punched  Cards 

Card  #1  (TITLE  (K),  K  =  3,20) 


Format 

(18A4) 


TITLE  —  72  column  title  for  plot.  (The  program  defines  the  first  word  of  the  title  as  either 
SURFACE  or  BOTTOM.) 


1.  Punched  Cards 
Card  #2  KTAPE,  MRG 
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Format 

(215) 


KTAPE—  Logical  unit  number  of  input  file. 

MRG  —  maximum  range  (nmi)  for  which  transmission  loss 
is  to  be  computed. 

Card  #3  XM1N,  XMAX,  YMIN,  YMAX,  XAX1S,  YAXIS  (6F10.5) 


XMIN  - 
XMAX  - 
YMIN  - 
YMAX  - 
XAXIS  - 
YAXIS  - 


minimum  x-coordinate  (initial  range  in  nmi). 

maximum  x-coordinate  (final  range  in  nmi). 

bottom  y  coordinate  (maximum  transmission  loss  in  dB). 

top  y  coordinate  (minimum  transmission  loss  in  dB). 

length  of  x-axis  in  inches. 

length  of  y-axis  in  inches  (<  10  in.). 


2.  File:  Assign  output  file  from  CONTOUR  to  FTONFOOl  where  N  =  KTAPE  (Card  #2). 
Subroutines  Called:  PLOTAX,  PLOTPT 


Output 


1.  Printed: 

Transmission  loss  in  dB  units  as  a  function  of  horizontal  range  in  steps  of  0.5  nmi  is  listed 
in  a  two  column  table. 


2.  Plot: 

A  Calcomp  plot  of  transmission  loss  vs  range  is  generated. 

Program  Name:  BISTRV 

Purpose:  Calculate  either  monostatic  or  bistatic  reverberation  as  a  function  of  time.  Compute  the  ver¬ 
tical  distribution  of  reverberation  at  up  to  20  specified  times.  If  the  source  and  receiver  are 
separated  in  horizontal  range  (bistatic  case)  then  the  calculation  is  correct  only  if  sound 
velocity  and  bottom  are  range  independent.  Either  surface  or  bottom  reverberation  is  calcu¬ 
lated  in  a  single  execution  of  program  BISTRV. 


Input 

1 .  Punched  Cards 

Format 

Card  #1  TITLE 

TITLE  —  80  column  title  for  printed  output. 

(20A4) 

Card  #2  MONO,  MXSORD,  MXRORD 

MONO  —  flag  for  type  of  calculation. 

(315) 

—  0  if  a  bistatic  calculation  is  to  be  made, 

=  1  if  a  monostatic  calculation  is  to  be  made. 


MXSORD  —  maximum  order  of  contour  to  use  from  source  file. 

MXRORD  —  maximum  order  of  contour  to  use  from  receiver  file. 


25 


1.  Punched  Cards 
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Format 

Card  #3  RABSMA,  PABSMA  (2F10.2) 

RABSMA  —  maximum  range  from  receiver  (nmi)  to  consider  in 
reverberation  calculation. 

PABSMA  —  maximum  range  from  source  (nmi)  to  consider  in 
reverberation  calculation. 

Card  #4  X1SOUR,  X1REC,  DXPR,  SEP  (4F10.2) 

X1SOUR  —  minimum  p  coordinate  (source)  of  grid  point  (nmi). 

X1REC  —  minimum  r  coordinate  (receiver)  of  grid  point  (nmi). 

DXPR  —  range  increment  for  p  and  r  coordinates. 

SEP  —  source-receiver  separation  in  horizontal  range  (nmi). 


Note  1:  For  monostatic  calculation,  input  variables  XIREC  and  SEP  are  disregarded.  They  are 
defined  by  the  program  as  explained  below. 


Note  2:  The  variables  defined  on  Card  #4  refer  to  a  ( p,r )  coordinate  system  where  p  is  the  hor¬ 
izontal  range  from  the  source  S  to  the  scattering  point  P  and  r  is  the  horizontal  range 
from  P  to  the  receiver  point  R.  The  distances  p,r  and  SEP  are  shown  in  Fig.  9.  Here 
the  x-j  plane  is  the  plane_ containing  the  ocean  boundary  (surface  or  bottom),  Fis  the 
scatter  point,  and  S  and  R  are  projections  of  S  and  R  onto  the  boundary  plane.  It  is 
apparent  that  for  a  given  S-R  separation  of  SEP,  the  following  three  inequalities  must 
be  satisfied: 


p  +  r  ^  SEP, 
r  +  SEP  >  p, 


p  +  SEP  >  r. 


v 


Fig.  9  —  Notation  and  geometry  of  bistatic  configuration 
for  one  scattering  path 
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Figure  10  shows  the  scattering  region  in  the  p-r  plane  along  with  other  variables 
defined  on  Card  #4. 

There  are  values  of  certain  parameters  that  will  collapse  the  general  reverberation  calcu¬ 
lation  to  the  monostatic  case.  These  are: 

1.  X1SOUR  =  X1REC  (usually  0) 

2.  SEP  =  DXPR/2 


Fig.  10  —  Bistatic  scattering  region  in  p-r  plane 


Under  these  conditions  the  ( p ,  r)  grid  points  are  such  that  p  =  r  as  shown  in  Fig.  1 1 . 
The  program  automatically  defines  X1REC  equal  to  XI  SOUR  and  SEP  equal  to 
DXPR/2  for  monstatic  calculation  (MONO  =  1). 


Fig  11  —  Monostdiii.  suiMering  region  in  />-f  plane 


i  - 
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There  are  requirements  which  must  be  observed  when  calculating  bistatic  reverberation. 
For  a  given  source  receiver  separation  SEP,  one  must  satisfy: 

1.  DXPR  “  SEP/n,  where  n  is  a  positive  integer, 

2.  tXISOUR  -  XIRECi  =  (m  +  1/2)  (DXPR),  where  m  is  a  nonnegative 
integer. 

To  illustrate,  the  actual  selection  of  parameters  for  a  bistatic  calculation  would  probably 
occur  in  the  following  way: 

Step  1:  Select  grid  increment  DXPR  based  upon  the  given  separation 
SEP  such  that  DXPR  —  SEP/n. 


Xj?\ 


Step  2:  Select  X1SOUR  -  0,  X1REC  =-  DXPR  /  2  (or  vice  versa). 
The  Jacobian  of  the  transformation  between  x—y  coordinates 
and  p—r  coordinates  is  not  defined  on  the  boundaries 
of  the  scattering  region  (see  Fig.  10).  A  sufficient 
condition  for  grid  points  to  miss  the  boundaries  is  that 
X1REC  not  be  an  integer  multiple  of  DXPR.  The  further 
restriction  that  X1REC  =  DXPR/2  assures  that  the  grid 
squares  represent  an  area  equal  to  the  area  of  the  integration 
region  (except  for  very  early  arrivals). 

For  example  choosing  DXPR  —  SEP/4  (Step  1)  along  with  X1REC  ■ 
**  DXPR/2  (Step  2)  yields  the  grid  shown  in  Fig.  10. 


0  and  XI SOUR 


Steps  1  and  2  determine  the  grid  to  calculate  reverberation  for  all  times.  If  reverbera¬ 
tion  is  desired  only  for  scattering  areas  beyond  certain  ranges  or  for  times  greater  than  a 
minimum  time,  then  X1SOUR  and  X1REC  may  be  redefined  as  follows.  In  Step  2, 
select  X1SOUR  =  mi  (DXPR),  X1REC  =  (m2  +  1/2MDXPR)  (or  vice  versa)  where 
mi  and  m2  are  integers  chosen  so  that  XI SOUR  and  X1REC  are  equal  to  the  minimum 
ranges  to  be  considered  for  the  reverberation  calculation. 

If  reverberation  is  desired  for  times  greater  than  some  minimum  time  t0,  then  appropri¬ 
ate  values  of  X1SOUR  and  X1REC,  can  be  defined  by  taking  the  velocity  of  sound  in 
the  ocean  to  be  approximately  0.8  nmi/s.  Remember  that  the  time  at  which  reverbera¬ 
tion  is  calculated  is  the  two-way  travel  time  from  the  source  to  the  scattering  point  and 
then  to  the  receiver.  To  illustrate  this  last  point,  suppose  reverberation  values  are 
desired  for  times  greater  than  90  s,  which  corresponds  to  a  total  distance  of  approxi¬ 
mately  72  nmi.  One  would  then  choose  mi  -  m2  such  that  XI SOUR  and  X1REC  are 
somewhat  less  than  36  nmi  (say  30  nmi),  thus  eliminating  calculations  for  reverbera¬ 
tion  times  not  desired.  A  similar  calculation  would  apply  to  the  monostatic  case. 


Card  #5  DUR 

DUR  —  pulse  duration  of  source  (s). 


(F10  7) 


$ 
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1.  Punched  Cards  Format 

Card  #6  Tl,  DT,  T2,  THETA1,  DTHETA,  THETA2  (6F10.2) 


T1 

DT 

T2 

THETA1  - 

DTHETA  — 
THETA2  _ 


initial  time  at  which  to  compute  reverberation  (s) 
time  increment  (s) 

final  time  at  which  to  compute  reverberation  (s) 

initial  vertical  angle  (degree)  at  which  to  compute  the  vertical  distribution  of 
reverberation  (-90°  direction  is  toward  surface), 
angle  increment  (degree). 

final  vertical  angle  (degree)  at  which  to  compute  the  vertical  distribution  of 
reverberation. 


Note  1:  A  maximum  of  2000  times  may  be  defined. 

A  maximum  of  181  angles  may  be  defined. 

If  the  vertical  distribution  of  reverberation  is  not  to  be  calculated,  the  last 
three  fields  in  this  card  may  be  left  blank 

Note  2:  The  reverberation  at  time  t  is  defined  to  include  all  contributions  to  rever¬ 
beration  occurring  in  the  time  interval  from  t  to  t  +  DT.  This  value  is  then  divided 
by  DT  so  that  the  final  value  is  actually  a  reverberation  density  per  unit  time. 

Type  #7  VERTTI(IV)  (8F10.2) 


VERTTI— the  IVth  time  at  which  to  compute  the  vertical  distribution 
of  reverberation.  Up  to  20  times  are  allowed  with  a  zero 
value  indicating  the  end  of  input  times. 

Note:  At  least  one  card  is  required  (blank  card  indicates  no  vertical  distribution  calculation). 
The  reverberation  contribution  at  angle  0  at  time  t  is  determined  by  that  portion  of  the 
the  total  contribution  to  reverberation  density  (defined  above)  corresponding  to  rays 
having  receiver  vertical  angles  in  the  angle  interval  centered  at  0  and  having  angular 
size  DTHETA.  The  value  obtained  is  divided  by  DTHETA  to  yield  reverberation  den¬ 
sity  per  unit  time  per  unit  angle. 

Type  #7  SIG  (8F10.2) 

SIG  —  An  array  of  91  values  defining  a  backscattering  table  for 
the  appropriate  boundary  as  a  function  of  grazing  angle 
from  0°  (tangent  to  boundary)  to  90°  (normal  to  boundary) 
in  1°  steps.  Values  are  in  dB  <  0. 

2.  Files:  Assign  source  file  generated  by  program  CONTUR  TO  FTO1F001. 

Assign  receiver  file  generated  by  program  CONTUR  to  FT02F001. 


FRANCHI,  GRIFFIN,  AND  KING 


Subroutines  Called:  There  are  no  subroutines  called. 


Output 

1.  Printed. 

a.  Listing  of  input  parameters. 

b.  Table  of  reverberation  in  both  dB  units  and  intensity  units 
vs  time.  The  reverberation  values  are 

relative  to  an  on  axis  source  level  of  0  dB. 

c.  Tables  of  the  vertical  distribution  of  reverberation  at  the  receiver 
vs  angle  for  the  times  specified  (if  calculated). 

2.  Files:  FT10F001  Catalogue  this  file  for  use  as  input  file  for 

program  AVREVB.  The  file  contains  the  reverberation  vs  time 
results  just  as  they  are  printed. 

FT20F001  This  file  contains  the  vertical  distributions  of  reverberation 

Of  calculated). 


Program  Name:  AVREVB 

Purpose:  Compute  a  weighted  average  of  one  or  more  reverberation  envelopes  calculated  by  program 
BISTRV.  The  resulting  curve  is  time  averaged,  adjusted  for  source  and  noise  levels,  and 
plotted. 

Usage: 


Format 

(20A4) 

(315) 


NOCRVS  —  number  of  envelopes  to  be  averaged. 
NOAVE  —  number  of  discrete  times  that  time  averaging 
is  performed  over. 

NOPTS  —  number  of  points  in  each  envelope. 


Input 

1.  Punched  Cards 

Card  #1  TITLE 

TITLE  —  80  column  title  for  plot. 


Card  #2 


NOCRVS,  NOAVE,  NOPTS 


Card  #3  XMIN,  XMAX,  YBOT,  YTOP,  XAXIS,  YAXIS,  ANL,  SL 


XMIN  —  minimum  time  (s)  used  for  labeling  and  scaling  plot. 

XMAX  —  maximum  time  (s)  used  for  labeling  and  scaling  plot. 

YBOT  —  minimum  reverberation  level  (dB)  used  for  plotting. 

YTOP  —  maximum  reverberation  level  (dB)  used  for  plotting. 

XAXIS  —  length  of  x  or  time  axis  in  inches. 

YAXIS  —  length  of  y  or  reverberation  level  axis  in  inches  10) 

ANL  —  noise  level  in  dB  re:  source  level. 

SL  —  source  level  in  dB  re:  1  yard 


(8F10.3) 
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1.  Punched  Cards 
Type  #4  WT 


Format 

(F5.0) 


WT  —  The  weight  assigned  to  the  Nth  envelope  input,  N  =  1,  2,  3, . NOCRVS, 

2.  Files.  Assign  output  files  from  B1STRV  to  FT01F00N  where  N  =  1,  2,  3  ....  NOCRVS.  Thus 
BISTRV  is  run  "NOCRVS"  times,  creating  "NOCRVS"  files  each  of  which  must  be  separately  assigned 
before  running  AVREVB. 

Subroutines  Called.  PLOTAX,  PLOTPT 

Output 

1.  Printed: 

a.  Input  variables  except  for  plotting  parameters. 

b.  Reverberation  values  in  dB  both  before  and  after  being 
adjusted  by  source  and  noise  levels  and  being  time 
averaged. 

2.  Plot:  A  Calcomp  plot  of  the  adjusted  reverberation  vs  time  is  generated. 

EXAMPLES  OF  REVERBERATION  CALCULATIONS 

Two  examples  illustrating  the  use  of  the  sequence  of  computer  programs  for  predicting  acoustic 
reverberation  are  presented.  The  first  example  is  a  "pseudo"  monostatic  reverberation  calculation,  i.e., 
a  calculation  where  the  acoustic  source  and  receiver  positions  differ  in  depth  but  not  in  horizontal 
range.  The  second  example  is  a  bistatic  reverberation  calculation. 

Figures  12  and  13  show  how  the  computer  programs  are  used  to  calculate  reverberation.  Input  for 
program  PROFIL  consists  solely  of  punched  cards.  Inputs  for  all  other  programs  consist  of  punched 
cards  and  data  files  created  by  previously  executed  programs. 


SURFACE 

REVERBERATION 


.  12  —  Implementation  of  ray-tracing  and  reverberation  programs 
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Fig,  13  —  Implementation  of  MERGE  program 


The  computer  programs  were  compiled  and  executed  on  the  Texas  Instruments  Advanced 
Scientific  Computer  (ASC)  at  the  U.  S.  Naval  Research  Laboratory. 

Example  I:  "Pseudo"  Monostatic  Reverberation 

Ocean  Environment 

Range  independent  sound  speed  field  (North  Atlantic  winter  profile  used). 

Linearly  upward  sloping  ocean  bottom  with  slope  of  lm/km. 

Bottom  loss  table  and  backscattering  tables  defined  consistent  with  area  of  North  A  .itic 
corresponding  to  sound-speed  profile 

Acoustic  Source  Characteristics 

Source  depth:  400  m. 

Source  frequency:  200  Hz. 

Source  beam  pattern:  omnidirectional  in  azimuthal  direction,  possesses  vertical  directivity  (see 
Fig.  15). 

Pulse  length:  10  s. 

Source  level:  200  dB. 

Acoustic  Receiver  Characteristics 
Receiver  depth:  500  m. 

Receiver  beam  pattern:  omnidirectional  both  horizontally  and  vertically. 

Calculations  Desired 

Acoustic  reverberation  from  both  the  ocean  surface  and  ocean  bottom  as  a  function  of  time  up  to 
400  s  after  onset  of  source  pulse. 

Figure  14  shows  the  output  plot  of  program  PLTENV.  Although  the  velocity  profile  is  not  shown 
extending  to  the  line  representing  the  ocean  bottom,  the  profile  is  linearly  extrapolated  to  the  bottom 
for  all  ray  tracing  calculations. 
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Fig.  14  —  Output  of  PLTENV  (Example  1) 


Ray  tracing  calculations  were  performed  by  program  REVRAP.  The  resulting  rays  were  reordered 
in  programs  PLTCON  and  CONTUR  to  produce  order  contours.  Surface  and  bottom  order  contours 
were  formed  and  plotted  by  program  PLTCON.  Order  contours  for  the  source  depth  are  shown  in  Figs. 
15  and  16.  Similar  plots  would  be  produced  for  the  receiver  depth. 

The  angle  fans  used  in  program  REVRAP  for  the  source  and  receiver  were  motivated  by  noting 
that  the  greatest  variations  in  intensity  at  the  boundaries  and  the  major  contributions  to  reverberation 
from  the  boundaries  will  occur  for  rays  whose  initial  angles  lie  in  an  interval  slightly  larger  than  the 
interval  defined  by  the  surface  grazing  ray  angle  and  the  bottom  grazing  ray  angle.  Using  Snell’s  law 
and  extrapolating  the  input  velocity  profile  (or  the  earth’s  curvature  corrected  profile  printed  by  pro¬ 
gram  PROFIL)  to  the  bottom,  the  bottom  grazing  angles  and  the  surface  grazing  angles  were  approxi¬ 
mated.  The  actual  values  calculated  for  the  source  depth  were  an  initial  angle  of  4.8°  for  grazing  to  the 
surface  and  initial  angles  of  13.7°  at  a  bottom  depth  of  4700  m  and  12.8°  at  a  bottom  depth  of  4329.6  m 
(the  depth  at  maximum  range)  for  grazing  to  the  bottom.  From  these  estimates,  the  angles  chosen  for 
program  REVRAP  were:  1.0°  steps  from  -45°  to  -20°,  0.5°  steps  from  -19.5°  to  -15.0°,  0.2°  steps 
from  -14.8°  to  0°,  and  similar  increments  for  the  positive  angles. 


Rays  were  traced  to  a  maximum  range  of  200  nmi.  If  one  assumes  that  acoustic  waves  in  water 
travel  at  a  rate  of  approximately  0.8  nmi/s  then  200  nmi  corresponds  to  a  travel  time  of  250  s.  The 
maximum  two-way  travel  time  from  source  to  scattering  point  and  back  to  receiver  could,  therefore,  be 
as  large  as  500  s,  which  is  well  above  the  400  s  defined  by  the  example. 

Orders  up  to  order  18  were  requested  in  program  REVRAP  to  ensure  that  all  significant  rays 
needed  for  intensity  calculations  were  traced.  The  maximum  reverberation  time  of  400  s  corresponds 
to  a  source/receiver  to  scattering  point  range  of  about  160  nmi.  Examination  of  Fig.  15  shows  that  the 
choice  of  18  as  a  maximum  order  (17  for  surface  orders)  was  adequate. 
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Fig.  16  —  Order  contours  for  source  to  bottom  (Example  I) 


The  vertical  beam  pattern  defined  for  the  source  is  shown  in  Fig.  17.  This  pattern  was  input  in 
program  CONTUR  by  means  of  punched  cards.  Note  that,  since  the  ray  angle  fans  in  program  RF.V- 
RAP  were  defined  from  —45°  to  45°,  values  of  the  beam  pattern  outside  that  interval  were  not  used  in 
any  computations;  i.e.,  the  angular  limits  of  the  vertical  beam  pattern  are  determined  by  the  ray  angle 
fan.  Program  CONTUR,  which  calculates  final  ray  intensities,  was  executed  four  times,  twice  with  the 
beam  pattern  defined  for  the  source  and  twice  with  the  omnidirectional  pattern  defined  for  the  receiver 
(see  Fig.  ,12)  The  two  executions  for  each  beam  pattern  produce  first  surface  and  then  bottom  order 
contours. 


NRL  REPORT  8721 


Fig.  17  —  Source  vertical  beam  pattern  (both  examples) 


Results  of  four  executions  of  program  TOTLOS  are  shown  in  Figs.  18  through  21.  The  greater 
transmission  loss  for  the  source  (400-m  depth)  is  due  to  its  vertical  beam  pattern,  the  receiver  (500-m 
depth)  being  omnidirectional  in  the  vertical. 

Calculation  of  reverberation  was  accomplished  by  program  B1STRV.  Two  executions  were 
required,  one  for  surface  re  verbena  ion  and  one  for  bottom  reverberation.  Program  AVREVB  added 
the  source  level  of  200  dB  to  each  reverberation  calculation  and  created  a  plot  of  reverberation  vs  lime. 

Surface  and  bottom  reverberation  curves  are  shown  in  Figs.  22  and  23.  It  is  expected  that  regions 
of  high  reverberation  generally  correspond  to  regions  of  low  transmission  loss.  This  is  indeed  the  case 
as  can  be  observed  by  converting  the  reverberation  arrival  time  axes  of  Figs.  22  and  23  to  range  from 
source/receiver,  and  then  comparing  Figs.  22  and  23  with  Figs.  18  through  21. 
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Example  2:  Bistatic  Reverberation 

This  example  differs  from  Example  1  in  that  the  ocean  bottom  is  now  defined  as  flat  and  the 
receiver  is  separated  from  the  source  horizontally  by  40  nmi,  thus  defining  a  true  bistatic  situation. 

Input  parameters  used  for  all  computer  program  executions  with  the  exception  of  programs  PRO- 
FIL  and  BISTRV  were  exactly  the  same  as  those  used  in  Example  1.  The  ocean  bottom  was  defined  to 
be  of  constant  depth  in  program  PROFIL  and  the  40-nmi  source-receiver  separation  was  input  in  pro¬ 
gram  BISTRV. 

Output  plots  of  program  AVREVB  are  shown  in  Figs.  24  and  25.  The  slight  rippling  of  the  curves 
is  due  to  a  grid  size  of  1.0  nmi  being  defined  in  BISTRV  (variable  DXPR)  rather  than  a  grid  size  of  0.2 
nmi  which  was  used  in  Example  1.  The  smaller  grid  size  would  have  eliminated  the  ripples.  However, 
the  associated  execution  time  for  BISTRV  would  have  increased  more  than  tenfold  over  the  actual  exe¬ 
cution  time  using  the  larger  grid  size. 

Execution  Time 

The  solution  to  each  of  the  two  examples  required  18  separate  program  executions  as  shown  in 
Fig.  12.  The  total  time  required  to  execute  all  necessary  programs  for  Example  1  was  approximately 
200  s.  A  time  of  125  s  was  required  for  the  two  executions  of  program  REVRAP  alone.  Example  2 
took  about  570  s  to  run,  410  s  of  which  were  required  to  run  program  BISTRV  twice.  Table  2  shows 
typical  execution  times  for  each  program. 


Bistatic  reverberation  envelope  for  surface  interaction 
(Example  2) 


Fig.  25  —  Bistatic  reverberation  envelope  for  bottom  interaction 
(Example  2) 


Table  2  —  Execution  Time  for  a  Single  Execution 
of  Computer  Programs  used  for  Examples  1  and  2 
using  the  Texas  Instruments  ASC  Computer  at  NRL 


Program  Name 

Typical  Execution  Time  (s) 

PROFIL 

0.2 

PLTENV 

0.5 

revrap 

62 

MERGE 

_ 

PLTCON 

5 

CONTUR 

6 

TOTLOS 

1 

BISTRV  (Example  2) 

Surface 

250 

Bottom 

160 

BISTRV  (Example  1) 

Surface 

23 

Bottom 

15 

AVREVB 

1 
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Appendix 

TRANSFORMATION  TO  BIRADIAL  COORDINATES 


We  wish  to  transform  Eq.  (15), 

£  I  Cv„ur.n 

,imt.X.Y)  jtniX.Y) 


B(0„  <*>,)  BU>,,  $,) 
L,  (X,Y)  Lj(X,  Y) 


«r(#;,  ■*>,, *,)  t/y.  (ad 


into  the  biradial  coordinate  system.  The  individual  symbols  in  Eq.  (Al)  are  defined  in  the  main  text. 
For  brevity  we  write 


RU)  =  /,  n  f  dX  dY. 


Since  we  have  assumed  a  flat  bottom  and  a  sing'^  sound-speed  profile  for  the  bistatic  case,  propa¬ 
gation  characteristics  are  symmetric  about  the  vertical  plane  containing  the  source  and  receiver.  Also, 
in  this  report  both  the  source  and  the  receiver  beam  patterns  are  omnidirectional  in  the  horizontal 
plane.  Therefore  /  is  symmetric  about  the  line  joining  the  (X,  Y)  coordinates  of  the  source  and 
receiver.  This  line  is  referred  to  as  the  source-receiver  axis.  Thus, 


//,  f  dx  dy  =  2  / SA+fdx  dy- 


where  A  +  is  the  upper  half  plane  of  the  boundary  area  A  about  the  source-receiver  axis.  Clearly,  a 
more  general  formulation  where  the  beam  patterns  are  not  symmetric  is  possible,  but  for  simplicity,  we 
restrict  our  attention  to  the  symmetric  case. 

If  we  refer  back  to  Fig.  8,  we  can  see  that  the  Cartesian  coordinates  and  the  biradial  coordinates 
are  related  by  the  equations 

p2  =  Y2  +  {h  +  X)2  ( A4) 


r2  -  Y2  +  (h  -  X)2. 


Solving  for  X(p,  r)  and  Y(p,  r ),  we  obtain  the  transformation  equations 


[8 h2(r2  +  p2)  -  (r2- p2)  -  I6h 


From  Eqs.  (A4)  and  (A5)  it  is  clear  that  each  point  in  A  has  a  unique  image  point  in  the  (p,  r)  plane. 
What  we  seek  now  is  the  image  of  A+  in  the  (p,  r)  plane.  This  is  the  set  of  ( p ,  r)  such  that 

[8 h2(r2  +  p2)  -  (r2  -  p2)2  -  16 h4)  >  0.  (A8) 


NRL  REPORT  8721 


Factoring  inequality  (A8)  we  obtain 

(r  —  p  +  2h)  (r  +  p  +  2 h)  (p  —  r  +  2 h)  (r  +  p  —  2 h)  >  0. 


(A9) 


An  examination  of  this  inequality  shows  the  set  of  (p,  r)  we  seek  is  the  region  K  shaded  in  Fig.  8,  and 
that  the  source-receiver  axis  Y  -  0  folds  into  the  boundary  of  K. 


The  transformation  of  Eq.  (A2)  into  the  biradial  coordinate  system  is  thus 

/  L fdx  dy  =  2f  L+ fdx  dy  =  2  /  SK  /|y|</r  dp- 

We  now  turn  our  attention  to  the  evaluation  of  the  Jacobian  J  of  the  transformation 

bX_  BY 
dr  dr 
dX  d  Y 

dp  dp 

Substituting  the  appropriate  partial  derivatives  into  (All)  yields 


(A10) 


(All) 


—r 

r(4h}—  r2  +  p2) 

2h 

8  h2Y 

J_ 

p(4h2  +  r2  —  p2) 

2  h 

ih2Y 

(A12) 


where  Fhas  been  left  in  the  partial  derivatives  for  brevity.  Solving  for  J  and  simplifying  we  obtain 


Ul--^ 


2si_ 


4 hY  (8 hHr2  +  p2)  -  (r1  —  p2)2  -  16 A4! 1/2 


(A13) 


Clearly  J  is  nonvanishing  and  well-defined  except  where  Y  -  0.  As  was  shown  in  the  derivation  of  the 
region  of  integration  K,  this  corresponds  to  the  boundary  of  K.  Thus  the  transformation  is  well-defined 
inside  of  K  and  Eq.  (Al)  transforms  to 


c  c  <}> i)  3  (,0ii  0/) 

R ’  2,1  f  J*  ,<&,  C‘-uM  r) 


(A  14) 


